﻿/*************************
用CRS法表示稀疏矩阵
CPP文件
*************************/
#include<iostream>
#include "lecture6.h"
using namespace std;
//const int NN=1000;

void SparseMatrix::init()
{
	cout<<"行数:";cin>>row;cout<<endl;
	cout<<"列数:";cin>>col;cout<<endl;
	cout<<"行数=="<<row<<'\t'<<"列数=="<<col<<endl;
	cout<<"==================四阶矩阵的LU分解求解========="<<endl;
	cout<<"输入矩阵:"<<endl;
	int i;
	for(i=0;i<row;i++)
	{
		for(int j=0;j<col;j++)
		{
			cin>>spm[i][j];
		}
	}
	cout<<"输入Y值:"<<endl;
	for(i=0;i<row;i++){
		cin>>yy[i];
	}
}

void SparseMatrix::compressmatrix()
{
	int count=0;
	int fcount=0;
	int i,j;
	for(i=0;i<row;i++){
		for(j=0;j<col;j++){
			if(spm[i][j]==0){
				continue;
			}
			else
			{
				val[count]=spm[i][j];
				colidx[count]=j;
				rowptr[count]=count;
				count++;
			}
		}
	}	
}

void SparseMatrix::crsdisplay()
{
	int count;
	int i;
	for(count=0;count<NN;count++){
		if(val[count]==0){
			break;
		}
	}
	for(i=0;i<count;i++){
		crsm[0][i]=val[i];
		crsm[1][i]=colidx[i];
		crsm[2][i]=rowptr[i];
	}
	cout<<"CRS矩阵:"<<endl;
	for(i=0;i<3;i++){
		for(int j=0;j<count;j++){
			cout<<crsm[i][j]<<'\t';
		}
		cout<<endl;
	}
}

void SparseMatrix::LUsolve()
{
	float l[NN][NN]={0};
    float u[NN][NN]={0};
    float y[NN]={0}; //定义数组Y
    float x[NN]={0}; //定义数组X
    float sum=0;
    int i,j,k;
    /*初始化矩阵l*/
    for(i=0; i<NN; i++)
    {
        for(j=0; j<NN; j++)
        {
            if(i==j) l[i][j] = 1;
        }
    }
    /*开始LU分解*/
    /*第一步：对矩阵U的首行进行计算*/
    for(i=0; i<NN; i++)
    {
        u[0][i] = (float)(spm[0][i]/l[0][0]);
    }
    /*第二步：逐步进行LU分解*/
    for(i=0; i<NN-1; i++)
    {
        /*对“L列”进行计算*/
        for(j=i+1; j<NN; j++)
        {
            for(k=0,sum=0; k<col; k++)
            {
                if(k != i) sum += l[j][k]*u[k][i];
            }
            l[j][i] = (float)((spm[j][i]-sum)/u[i][i]);
        }
        /*对“U行”进行计算*/
        for(j=i+1; j<col; j++)
        {
            for(k=0,sum=0; k<row; k++)
            {
                if(k != i+1) sum += l[i+1][k]*u[k][j];
            }
            u[i+1][j] = (float)((spm[i+1][j]-sum));
        }
    }

    /*回代方式计算数组Y*/
    for(i=0; i<row; i++)
    {
        y[i] = gety(l,yy,y,i);
    }

    /*回代方式计算数组X*/
    for(i=row-1; i>=0; i--)
    {
        x[i] = getx(u,y,x,i);
    }

    /*显示数组X*/
    //printf("\n\n数组X：\n");
    cout<<"数组X："<<endl;
    for(i=0; i<row; i++)
    {
        //printf("x%d = %0.3f\n", i+1,x[i]);
        cout<<"x["<<i+1<<"]="<<x[i]<<endl;
    }
	cout<<"============L矩阵========"<<endl;
	for(i=0;i<NN;i++){
		for(j=0;j<NN;j++){
			cout<<l[i][j]<<'\t';
		}
		cout<<endl;
	}
	cout<<"============U矩阵========="<<endl;
	for(i=0;i<NN;i++){
		for(j=0;j<NN;j++){
			cout<<u[i][j]<<'\t';
		}
		cout<<endl;
	}
}

float SparseMatrix::getmx(float a[NN][NN],float x[NN],int i)
{
    float mx = 0;
    int r;
    for(r=i+1; r<NN; r++)
    {
        mx += a[i][r] * x[r];
    }
    return mx;
}

float SparseMatrix::getmy(float a[NN][NN],float y[NN],int i)
{
    float my = 0;
    int r;
    for(r=0; r<NN; r++)
    {
        if(i != r) my += a[i][r] * y[r];
    }
    return my;
}

float SparseMatrix::getx(float a[NN][NN],float b[NN],float x[NN],int i)
{
    float result;

    if(i==NN-1) //计算最后一个x的值
        result = (float)(b[i]/a[NN-1][NN-1]);
    else //计算其他x值（对于公式中的求和部分，需要调用getmx()函数）
        result = (float)((b[i]-getmx(a,x,i))/a[i][i]);

    return result;
}

float SparseMatrix::gety(float a[NN][NN],float b[NN],float y[NN],int i)
{
    float result;

    if(i==0) //计算第一个y的值
        result = float(b[i]/a[i][i]);
    else //计算其他y值（对于公式中的求和部分，需要调用getmy()函数）
        result = float((b[i]-getmy(a,y,i))/a[i][i]);

    return result;
}


int main()
{
	SparseMatrix mx;
	mx.init();
	mx.compressmatrix();
	mx.crsdisplay();
	mx.LUsolve();

	return 0;
}
